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Abstract 

We consider the dynamics and kinetic roughening of wetting fronts in the case 
of forced wetting driven by a constant mass flux into a 2D disordered medium. 
We employ a coarse-grained phase field model with local conservation of density, 
which has been developed earlier for spontaneous imbibition driven by a capillary 
forces. The forced flow creates interfaces that propagate at a constant average ve- 
locity. We first derive a linearized equation of motion for the interface fluctuations 
using projection methods. Prom this we extract a time-independent crossover 
length £ x , which separates two regimes of dissipative behavior and governs the 
kinetic roughening of the interfaces by giving an upper cutoff for the extent of the 
fluctuations. By numerically integrating the phase field model, we find that the 
interfaces are superrough with a roughness exponent of \ = 1.35 ±0.05, a growth 
exponent of (5 = 0.50 ± 0.02, and £; x ~ v~ x l 2 as a function of the velocity. These 
results are in good agreement with recent experiments on Hele-Shaw cells. We 
also make a direct numerical comparison between the solutions of the full phase 
field model and the corresponding linearized interface equation. Good agreement 
is found in spatial correlations, while the temporal correlations in the two models 
are somewhat different. 



1 Introduction 

The dynamics and roughening of moving interfaces in a disordered medium is a subject 
of intense interest in non-equilibrium statistical physics |Tj. Examples where such 
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processes are relevant include thin film deposition [2j, slow combustion fronts in paper 
[3J, fluid invasion in fractals |4| and porous media 0112 El, and wetting and propagation 
of contact lines [HI El UHj- The understanding of the underlying physics involved in 
interface roughening is crucial to the control and optimization of these processes, with 
immediately apparent technological importance. Progress in the theoretical study of 
interface dynamics has been made over the last two decades and a number of theories 
have been developed [H E| which agree with the experimental findings in some cases. 
Much of the theoretical work has been based on analyzing (spatially) local interface 
equations, such as the celebrated Kardar-Paris-Zhang (KPZ) equation [TT], where the 
physics is governed by a (nonlinear) partial differential equation which couples the 
interface locally with itself and the quenched randomness. However, in many cases 
such a description is not possible [H IT2]. 

A particularly important class of problems in the field of kinetic roughening where 
local theories cannot provide a complete description are those involving fluid invasion 
in porous media, which are often experimentally studied using Hele-Shaw cells |1 3| 
EU HH1 EH1 or even paper as the wetting medium [13 EH E3 EDI- The reason for this 
is that if the transport of liquid to the advancing wetting front from the reservoir is 
neglected as in local theories, slowing down of the front in spontaneous imbibition of 
water in paper cannot be explained by local theories unless the liquid conservation 
law is included in some artificial way. To properly describe the dynamics of wetting 
fronts in random medium is a challenging task, and there are several recent theoretical 
attempts to this end [13 Ell El- In particular, in Refs. [TJl El E3 EE! Dube et 
al. developed a phase-field model explicitly addressing the issue of liquid conservation 
in the wetting of a random medium. This is achieved by a generalized Cahn-Hilliard 
equation with suitable boundary condition which couples the system to the reservoir. A 
variant of the sharp interface projection method [22EII was used to analytically obtain 
a non-local interface equation for the case of spontaneous imbibition, and from it a new 
time dependent length scale governing the kinetic roughening £ x ~ was extracted. 
In Refs. [I2J EH] the kinetic roughening of 2D wetting fronts was also analyzed, and 
estimates for the corresponding scaling exponents were obtained. Furthermore, in a 
recent comprehensive review paper [2B] the case of forced wetting in 2D was briefly 
discussed, and in Ref. [211 forced wetting in a 3D model of paper was also considered. 

In the present work our aim is to carry out a comprehensive analysis of wetting 
fronts in 2D in the case of forced wetting. To this end, we use the phase-field model 
of Ref. [12] with boundary conditions corresponding to a constant mass flux at the 
reservoir boundary. This makes the wetting fronts move at a constant average velocity, 
in contrast to the Washburn law for spontaneous wetting. We first expand on the 
standard projection method to make it usable under constant flow boundary conditions, 
and derive the linearized interface equations corresponding to the forced case. Analysis 
of these equations reveals a crossover length scale £ x , which is time-independent in 
contrast to the spontaneous wetting case [12], and depends on the interface velocity as 
£ x oc v _1 / 2 . It separates two regimes of dissipative behavior and governs the kinetic 
roughening of the interfaces by giving an upper cutoff for the correlation length of the 
interface fluctuations. By numerically integrating the phase field model, we find that 
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the interfaces are superrough with a roughness exponent of x — 1-35 ± 0.05, and a 
growth exponent of (3 = 0.50 ± 0.02. These results are in good agreement with recent 
experiments on Hele-Shaw cells [H]. We also make a direct numerical comparison 
between the solutions of the full phase field model and the corresponding linearized 
interface equation. Good agreement is found in spatial correlations, while the temporal 
correlations in the two models are somewhat different. 

This paper proceeds as follows: In Section [2] we define the phase field model, and 
present results from the projection and linearization procedure. We present the result- 
ing interface equations, as well as discuss how the crossover length scale £ x emerges. 
In Section [HI we present our numerical results for the driven interfaces in disordered 
medium. We consider both spatial and temporal correlation functions, as well as com- 
pare the numerical results from the linearized interface equation to the phase field 
model. Finally, we present our discussion and conclusions in Section |3J The Ap- 
pendices A and B contain some technical details on the projection and linearization 
procedures leading to the interface equations. 

2 Model for Wetting 

2.1 Definition of the Phase Field Model 

The model describes the dynamics of a liquid invading a disordered medium at a 
coarse-grained level. A phase field is used to describe the "wet" and "dry" phases, 
with a free energy functional such that the dimensionless phase field obtains the values 
= +1 and = — 1 at the wet and dry phases, respectively. Since the phase field is an 
effective density field it is locally conserved. Energy cost for an interface is added by 
the standard coarse-grained gradient squared term. The interaction energy between 
the random medium and the invading liquid is represented by a quenched random field 
linearly coupled to the phase field. This leads to the free energy density |T2j 

^[0(x, *)] = i(W(x, t)) 2 + V(0(x, t)) - a(x)0(x, t), (1) 

where V has two minima at = — 1 and = +1, and a is the quenched random field. 
The standard Ginzburg-Landau form is chosen for V, i.e. V((f>) = — 2 /2 + 4 /4 [12], 
and the quenched field obeys the relations 

(a(x)> = (2) 
(a(x)a(x')) = (Aa)^(x-x'). (3) 

The case of positive (negative) (a(x)) corresponds to the liquid spontaneously wetting 
(dewetting) the medium for the case of no external driving [T2] . 

The equation of motion for the conserved phase field is given by the continuity 
equation <9 t = —V ■ j and Fick's law j = — V//, where fi = 5 J- '/5<p. The result is 

d t 0(x, t) = V 2 Mx, t) = V 2 [-0(x, t) + 3 (x, t) - V 2 0(x, t) - a(x)] , (4) 
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which is essentially the Cahn-Hilliard equation [SOj. Note that dimensionless units have 
been set such that the constant relating the phase field gradient to the free energy and 
the mobility in Fick's law are both unity. This fixes the choice of units in Eq. (J2J. 

In Refs. [ElEU the case of spontanous, capillary driven wetting was modeled with 
Eq. (jl]) using boundary conditions where the chemical potential was set to a constant 
at the liquid reservoir y = 0. In order to model the experimental setup of driving the 
liquid from the reservoir, we define our phase-field in the half-plane {x|y > 0}, and at 
the line y = impose the boundary condition V/i = — Fy, where y is the unit vector, 
and F is a constant (flux) parameter (see Fig. 1(a)). On the top end of the system 
we set (f>(y — > oo) = —1, and use periodic boundaries in the x direction. Physically 
this corresponds to driving the liquid via a constant mass flow, leading to an interface 
propagating at a constant average velocity [28j. 

The initial condition for the phase field is given by a step function at some height 
H(0) = H(t = 0), 0(x,t = 0) = 1 - 2S(y - H(0)). H(0) is also a parameter in our 
model, but with the gradient boundary condition here its value is irrelevant, in contrast 
to the case of spontaneous imbibition, where H(Q) defines the initial average velocity 
of the interface [T2] . 

2.2 Linearized interface equation 

The quenched disorder field a(x) will cause an initially flat interface to kinetically 
roughen when it propagates as the liquid invades the medium. A key step in under- 
standing the physics of this process is writing an equation of motion for the ID single- 
valued height variable H(x,t), defined conveniently by the condition <p(x, H(x, t)) = 
as shown in Fig. 1(a). It is not obvious a priori that this can be done, since the 
problem is inherently non-local due to the conservation law [TJ], but it has been shown 
in Refs. [12] that this can be done using methods discussed in Refs. [SHUSHES]. In 
Appendix lAl we will discuss some technical details of the derivation. The main idea is 
to use the relevant Green's function of the problem defined by 

V 2 G(x,y\x',y') = S(x - x) 5(y-y'), (5) 

with appropriate boundary conditions. This leads to the integro-differential equation 
of motion 

2 J dx' G{x,H{x,t)\x',H{x',t)) dH ^ t ^ =a(x,H(x,t)) + o-K + A\ y=H ( x ,t), (6) 

where k is the interface curvature, and A is the boundary term, which is non- vanishing 
for inhomogenous boundary conditions. Note that Eq. (jHJ) holds for any geometry given 
the appropriate Green's function. For the half plane with von Neumann boundary 
condition at y — 0, the 2D Green's function is obtained by image charge method as 

G(x, y\x', y>) = ±- In [((x - x' f + (y - y'f ){{x - x'f + (y + y'f)] ■ (7) 
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The next step towards an explicit interface equation is to linearize in fluctuations 
around the disorder-free system solution, H(x,t) = H (t) + h(x,t), and transform to 
Fourier space, where the equations of different modes of fluctuation are decoupled. The 
half plane Green's function defined in Eq. (jZJ) is not square integratable, however, and 
thus it does not have a Fourier representation. However, we have found that we can 
avoid this problem by considering a finite strip of width L where < x < L, and derive 
a linearized equation of motion for the discrete Fourier modes of the fluctuations. Then 
we can take the limit L — > oo to obtain the equation of motion in Fourier space for an 
interface in half-plane geometry, and the end result is well-defined. This procedure is 
exposed in some detail in Appendix FBI 

The linearization procedure, by construction, gives a separate equation of motion for 
the mean interface position Ho(t), which turns out to be coupled to all the fluctuation 
Fourier modes h k , while the evolutions of the fluctuation modes are independent of 
each other. The resulting equations are given by 

H = j ; (8) 

^(l + e-W) = \ k \(-H h k (l-e- 2 W H °)-ak 2 h h + Vk (t)), (9) 
from which we immediately obtain the expected result that 

Ho(t) = y. (10) 

This should be contrasted with the Washburn law H (t) oc t 1 ^ 2 in the case of sponta- 
neous wetting [12]. It is interesting to note that the functional form of Eq. (0) for the 
fluctuations of the interface is similar to the case of spontaneous wetting in Ref. [12], 
except for some sign changes in the terms. It can be shown that these changes are a 
direct consequence of the differences between the Green's functions in the two cases, 
and they lead to significant differences in the behavior of the fluctuations as will be 
discussed below. 

The most immediate aspect of these equations is that the fluctuation equation is 
non-local in real space. This is to be expected due to the conservation law [12] ■ The 
locality of the interface equation in Fourier space owes to the fact that it has been 
linearized. The interface configuration couples to the disorder in a fundamentally non- 
linear manner, a fact that is somewhat obscured by the superficially simple form of the 
disorder term rjk, which is defined as 

*<*>-/ "<*.'» (") 

The moments of r\ k averaged over disorder realizations thus couple to the interface con- 
figuration realizations, which in turn are defined by the disorder. This makes the anal- 
ysis of rjk a formidable task that must be solved self-consistently involving the interface 
equation. In more explicit terms, in expressions such as (at(x, H(x, t))a(x', H(x', t'))), 
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the angular brackets denote averages over different realizations of H, and results can 
only be obtained numerically. 

In Eq. (jnj there are two terms that dissipate fluctuations corresponding to different 
physical effects: the surface tension term ak 3 hk and the liquid transport (conservation 
law) \k\Ho (l — e~ 2 ^ Ho ) h k . The surface tension dominates when kH <C ak 3 , leading 
to a time independent crossover length scale between the two terms given by 




This is in striking contrast to the spontaneous case, where the corresponding crossover 
length is time dependent with £ x oc t 1//4 [12]. Moreover, in the dispersion relation of 
the fluctuations there is an additional crossover in the transport term obtained by 
comparing the length scales H (t) and k. Namely, for kH (t) 3> 1 the transport term 
is kHo, while for kH (t) <C 1 it is 2k 2 H H , leading to a crossover between u oc k 
and uj oc k 2 in the dispersion relation of the interface fluctuations. Regardless of the 
magnitude of kH (t), we will show through numerical studies that the crossover length 
scale £ x controls the kinetic roughening of the interface in analogy to the spontaneous 
imbibition case, in that it defines an upper cut-off for fluctuations that are increasing 
in time and correlated by the surface tension. 



3 Numerical analysis 

The interface fluctuations in the presence of quenched disorder were analyzed by nu- 
merical integrations of Eq. (4) with the appropriate boundary conditions shown in the 
schematic representation of Fig. In order to implement the von Neumann boundary 
condition at y — 0, Eq. (4) is modified as follows: 

d t 0(x, t) = V 2 [-0(x, t) + 3 (x, t) - V 2 0(x, t) - FH (t) - a(x)] , (13) 

where the addition of the term FH (t) ensures that the interface will propagate at 
constant velocity v = H — F/2 due to the boundary condition d y fi\ y=0 = —F. The 
above equation is then solved with the Dirichlet boundary condition n\ y =o = 0. This is 
a numerical trick invoked to obtain a chemical potential profile consistent with the von 
Neumann boundary condition in the domain from the reservoir at y = to the interface 
at y = H. The position of the interface H(x,t) at each x was defined as H(x,t) = 
by linear interpolation between the points of the numerical grid. Without any loss of 
generality, the chemical potential at the boundary is chosen such that fi(x,y = 0) = 0, 
leading to 4>(x, y = 0) = 4>o, where <fio is the solution of — (fio + 4>l = FH (t) 1 . 

x It should be pointed out that the von Neumann boundary condition can also be implemented 
directly by using fj,(y = 0) = n{y = Ay) + FAy, where Ay is the size of the spatial discretization. 
The numerical implementation of the von Neumann boundary condition is then used to fix 4>(x, y = 
0) = (j)o, where 4>q is the solution of —fa + 0§ = ~ 'M^; V — Ay) + <p(x, y — Ay) 3 + FAy. We have 
compared the two different implementations and found no distinguishable differences. 
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Typical plots for the chemical potential and the density field along the y axis ob- 
tained from numerical integration of Eq. (jT3| without quenched disorder (Aa = 0) are 
shown in Figs. HJ(b) and (c), respectively. A set of successive interface configurations 
with v = 0.02, and Aa = 0.2 are also shown in Fig. [T] (d). It is interesting to note 
from Fig. QJc) that behind the moving interface the density field has a finite slope, 
which remains constant in time. This is due to the underlying local conservation law 
for (f>(x,y). We note that the slighly growing value of which is greater than +1 is 
just a numerical artifact and can be removed by using a method of moving box to hold 
the interface at a constant height, i.e. pulling down the disorder field at constant time 
intervals. 

3.1 Spatial roughness 

With Aa > 0, the driven wetting front kinetically roughens as can be seen in Fig. 
1(d). To characterize the spatial extent of the roughness, we first consider the spatial 
two-point correlation function 

G 2 (r,t) = ({h(x + r,t)-h{x,t)Y) l '\ (14) 

which is directly related to the structure factor S(k,t) = (hk(t)h-k(t)) . In the above 
equations the brackets denote an average over different configurations of random noise, 
and the overbar a spatial average over the system. 

In Fig. |2](a) we show numerical data for the spatial correlation function. We find 
that the correlation length of the roughness of the interface saturates after an initial 
growth. According to Eq. (|T2j) . the crossover length £ x is related to the interface 
velocity by £ x ~ v^ 1 ^ 2 . In the inset of Fig. [2(a) we plot the velocity dependence of the 
corresponding crossover length £ 2 found from G^r). We indeed find that this length 
£2 ~ f _0 ' 45 , which means that £2 oc £ x as in the case of spontaneous wetting, too [12] . 

An estimate for the global roughness exponent x can be obtained from the structure 
factor S(k,t), as shown in Fig. [21(b) . As expected, we find that S(k) ~ l/k 1+2x 
is well satisfied, with a global roughness exponent of x ~ 1-25, and a crossover to 
a plateau corresponding to distances larger than the intrinsic correlation length £ x , 
consistent with the analysis from the linearized interface equation. We actually found 
that the global roughness exponent slightly depends upon the velocity and increases 
with decreasing velocity until asymptotically approaching a value of about 1.35. For 
velocities v = 0.005, 0.002, and 0.001, with Aa = 0.1, it was found that x ~ 1-27, 1.35, 
and 1.37, respectively. Moreover, that global roughness exponent slightly depends on 
the strength of the noise. For example, for v = 0.005, x ~ 1-27 and 1.36 for Aa = 0.1 
and 0.2, respectively. We also compared the crossover lengths obtained from G 2 (r,t) 
to those from S(k,t), and found that the value of £ x from S(k,t) is about one-half of 
that from G 2 (r, t), independent of disorder strength. 

We also estimated the the local roughness exponent xioc according to the scaling 
relationship G 2 {r — l,t) ~ £ x ~ Xloc ~ u(xioc-x)/2 for different velocities [12], and found 
that Xioc ~ 1-0, as expected for superrough interfaces. The spatial correlation funnction 
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G%(r, t) should also follow the same scaling form G^O", t) ~ Aav x l 2 g{rv l l 2 ) as G^r, H) 
for spontaneous imbibition with fixed interface front height 1 1 2| . 

3.2 Temporal roughness 

To quantify the temporal development of the roughness, we consider the width of the 
interface defined by 



In the presence of quenched disorder the roughness initially increases as a power law of 
time, as shown in Fig. After a crossover, the roughness reaches a saturated regime. 
For small Aa = 0.1, only relatively low velocities were studied, because for velocity 
as high as v — 0.05, the roughness profile shows pronounced oscillations. The cause 
of such numerical oscillations was identified to be the numerical interpolation of the 
interface position between the grid points with time scale equal to lattice size over the 
velocity. It can be clearly seen from the inset of Fig. Elthat for the same noise strength, 
all the roughness curves follow the same initial growth profile, suggesting an universal 
growth exponent (3. It was found that for Aa = 0.1,0.2, and 0.3, the corresponding 
values of (3 are about 0.48, 0.50, and 0.52, respectively. 

The saturated width of the interface was found to be independent of the lateral 
system size as long as the lateral system size is larger than the intrinsic crossover length 
£ X j which has been derived from the linearized interface equation in the preceding 
section. It should be pointed out that in Ref. [2B|, where the case of driven wetting 
was briefly discussed, it is claimed that after an initial growth, the interface roughness 
follows a weak logarithmic growth. From Fig. El one can see that we do not find any 
evidence of such a logarithmic growth regime, although it could be too slow to be 
detected numerically. 

Assuming that the crossover length £ x controls the roughening process, we can use 
the results in Ref. |T2] and write a Family- Vicsek type of scaling relation 



Data collapse using this scaling from is presented in Fig. El Using the data collapse, 
we give our best estimates for the roughness and growth exponents as x — 1-35 ± 0.05 
and (3 = 0.50 ± 0.02. We note that Ref. [2Hj estimates that x ~ 1-25 and (3 « 0.4, 
with (3 — %/3. However, our results do not support this relation. 

In the case of spontaneous imbibition, it was found that the rough interfaces obey 
temporal multiscaling, with different scaling exponents for different moments of the 
time-dependent correlation functions ^2]. These functions and the corresponding ex- 
ponents are defined by 



w 2 (t) = ((h(x,t)-h(x,t)) 2 ). 



(15) 




(16) 



C q {t) = ([H(x, t + s)- H (t + s)- H(x, s) + Hois)}^, 



(17) 
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for q = 1, 2, 3, .... The correlation functions C^it) are shown in Fig. |IJa) for different 
velocities. The crossover time t x between the power law regime of C 2 {t) and saturation 
increases with decreasing velocities. The different functions C q {t) for q = 2,4 and 6 
are shown in Fig. Stb). At early times all the functions follow power law behavior 
C q (t) ~ t l3q , with exponents j3 q ~ 0.94 which are independent of q. This behavior 
is, however, observed only in time scales smaller than the disorder persistence time 
td = Ay/v, where Ay = 1 is the dimensionless spatial discretization step. If we 
consider C q (t) for t > td, we find evidence multiscaling with fi 2 ~ 0.79, 04 ~ 0.69 and 
/3q 0.54 in a small time regime, similar to the spontaneous case [12]. However, it is 
difficult to verify true multiscaling here because of the rapid crossover to the saturated 
regime, although we do expect avalanche type of motion to be present here, which 
often leads to multiscaling behavior |83] . 

3.3 Numerical Results from the Linearized Interface Equations 

An interesting question concerns the range of validity of the linearized interface equa- 
tion (LIE), Eq. (jHj, for the fluctuations of the height, in particular in the driven, 
nonlinear regime with noise included. This issue is also related to the possible exis- 
tence of universality classes of roughening for conserved systems with quenched noise; 
a problem for which there are virtually no analytical results. To this end, we have 
integrated Eq. (j^} numerically in time. To incorporate the nonlinear nature of the 
disorder T](x, H (x, t)), a back-and-forth Fourier transformation scheme is required at 
every time step. A square lattice landscape of independently Gaussian distributed 
noise was used for r)(x,y). We note that solving Eq. (jOJ) instead of the full 2D phase 
field model is numerically much easier. 

To compare the results with the phase-field model, we numerically computed the 
same set of correlation functions. From the analytic derivation of the interface equation, 
we can obtain a quantitative map between the parameters of the phase field model and 
the interface model, which is as follows. The interface velocities should obviously be 
the same. The disorder fields between the models are related by r] = Ma, where M 
is the mobility in Fick's law. In our dimensionless units M = 1. The effective surface 
tension in the phase field model is given by the standard form of the potential V(<f>) 
as a p f = J du (d u (f> (u)) 2 , where O is the ID kink solution of the disorder-free system. 
The surface tension in the interface equation comes out as a = Ma p f. 

The results from the LIE are collected in the insets of the corresponding phase field 
results in Figs. 2-4. Obviously, the length scale £ x is present in an identical manner in 
both cases. Remarkably enough, the corresponding results from the two cases are in 
most cases quantitatively close to each other. There are some important differences, 
however. First, the saturated interface widths differ by about 20%, the interfaces 
from the LIE being rougher. Since the interface model only takes into account linear 
dissipation effects, we would expect it to underestimate the stiffness of the interface, 
leading to larger roughness amplitudes than in the full phase field model. 

Next, we examine the data collapse for the interface roughness (Fig. 0J) using the 
scaling function of Eq. (fT6j) . The best collapse is obtained for slightly different values 
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of the exponents as compared to the phase-field model. The roughness exponent is 
close to the previous value, namely x — 1-27 ± 0.05. However, the growth exponent is 
now given by (3 = 0.37 ±0.04 for the LIE, where the difference to the phase field result 
is outside of the numerical uncertainties. In the LIE we also noted a slight deviation 
from the behavior of the crossover time t x oc as a function of the disorder strength 
Act. Another difference between the two models was found in temporal multiscaling. 
In the LIE each of the moments increases with a different exponent (3 q even at times 
smaller than the disorder persistence length td- The exponents are /?2 = 0.75 ± 0.03, 
/?4 = 0.55 ± 0.05, and {3q = 0.47 ± 0.03. The values of these exponents are lower than 
in the phase-field model, which is not surprising since the LIE has a lower value of (3, 
too. 

4 Summary and Conclusions 

In the present work, we have studied wetting of a disordered medium driven by a 
constant mass flux in a 2D system. Our model is a prototype phase field model 
incorporating mass conservation into the flow of two immiscible fluids, Eq. (|H). From 
this model we have derived non-local interface equations to lowest order in Fourier 
space fluctuations, Eqs. (jHJ) and 0. Because of the linearization, these equations are 
local in Fourier space. The constant flux boundary condition gives rise to a interface 
that moves with a constant velocity proportional to the flux. We have obtained a 
time-independent crossover length scale £ x oc \fo~jv from the interface equations. 
Numerically, we find that the kinetic roughening of an interface is governed by a scaling 
relation of the Family- Vicsek type, where £ x controls the extent of the fluctuations (for 
£ x < L) as given by Eq. (fT6|) . For the kinetic roughening of the interfaces, we find 
that they are superrough, with x = 1-35 ± 0.05 and xioc ~ 1- For temporal roughness, 
(3 = 0.50 ± 0.02, and there is numerical evidence of temporal multiscaling for t > td. 
We note that all these results are in qualitative but not in quantitative agreement with 
the case of spontaneous imbibition studied earlier in Refs. [12] for interfaces, which 
are kept at constant height in the steady-state regime described by Washburn's law. 

In addition to obtaining the LIE by analytic methods, we have also made a direct 
comparison between it and the full phase field model with quenched noise properly in- 
cluded. We find very good agreement between the spatial correlations of the interfaces, 
even including approximately the same roughness exponent of x ~ 1-3- However, the 
temporal correlations in the two cases are different: while the amplitude of the sat- 
urated roughness is larger in the LIE, but the growth exponent (3 = 0.37 ± 0.04 is 
smaller than the phase-field result (3 = 0.50 ± 0.02 in the phase field model. Also, 
spatial multiscaling is more clearly present in the LIE within our numerics. 

Our analytical and numerical results of forced wetting can be compared with ex- 
perimental results of kinetic roughening of an oil-air interface in a forced wetting where 
the experiments were done in a horizontal Hele-Shaw cell with quenched disorder |13j . 
It was found in the experiments that the growth exponent (3 ~ 0.5 which is nearly inde- 
pendent of the experimental parameters, and the roughness exponent x ~ 1-3, which, 
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however, depends on experimental parameters. While a fully quantitative comparison 
may be difficult, both exponents obtained experimentally are in very good agreement 
with our numerical results. Furthermore, the experiments confirm that the crossover 
length scales as the inverse of the square root of velocity, as found in our theory. Fur- 
ther experiments on the dependence of the results on other systems parameters would 
be most interesting. 
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Appendix 

A Projection to an Interface Equation 

In this Appendix we discuss the details of the sharp interface projection to obtain the 
interface equation, Eq. (jHJ) from the phase field equation of motion, Eq. (J3J). First we 
invert the phase field equation in a volume V with boundary S, by multiplying with 
the Green's function, integrating over V and applying Gauss's divergence theorem. 
The result is 

d 3 r'y/det(g')G(r, r')d t <j){r') = /i(r) + A, (18) 



v 



where the surface term A is explicitly given by 

A = J dS' • [G(r, r')W(r') - /i(r')VG(r, r% (19) 

The standard procedure is then to consider a single- valued sharp interface H (x, t), and 
transform to coordinates of distance along and perpendicular to this interface given by 
(s, u). In these coordinates the metric tensor is given by 
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1 

(1+uk) 2 



(20) 



where k is the curvature of the interface, defined via the unit tangent t and unit 
normal n of the interface as nt = <9 s n. The volume integration measure is the Jacobian 
J = ^det(g) = 1 + uk, and it must be positive definite. This limits the validity of the 
coordinates to the area not further from the interface than the radius of the interface 
curvature. 

Next, a number of standard approximations are made, including the small curvature 
approximation, which gives the Laplacian to first order in k as 

V 2^&_ ^ + (21) 

du 2 ds 2 du 
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For a sharp interface with small curvature, the phase field near the interface has the 
form given by the a = 0, ID kink solution 0o i n the normal direction, defined by 
dfyo — y'i^o)- The chemical potential is then fi ~ —Kd u (f)Q — a. Since the kink 
solution has a small gradient except at the interface, we can project Eq. (fTH| to 
the interface with the operator f dud u (j)o(u)[-}, and take the explicit sharp interface 
limit 0o — — 1 + 20(u). The projected equation involves contributions only from an 
area not further from the interface than the interface width, which is less than the 
interface radius of curvature by the virtue of the small curvature and sharp interface 
approximations. Therefore, the use of coordinates (s, u) is valid. Eq. ([THl) is then 
projected to 

2 J ds'G(s, 0\s', 0)d t u(s') = -ok - a(s, 0) + A| n=0 , (22) 

where a = \ J du (d u (j) (u)) 2 is the effective surface tension. This can be transformed 
to Cartesian coordinates using dsd t u = dxd t H(x,t), yielding Eq. (JHJ). 



B Linearization of the Interface Equation 

In this Appendix we describe in detail the method of using strip geometry to obtain the 
linearized interface equations (Eqns. JEJ) and JHJ)) from the full non-local sharp interface 
equation, Eq. Q. For the half-strip {(x,y)\x G [0,L],yE [0, oo]}, the Green's function 
for the Laplacian, with homogenous von Neumann boundary conditions at the strip 
edges, is given by 



G(x, y\x', y) = [\y - y'\ + y + y'\ 



1 1 /mrx\ 

> — cos — - — cos 

7r ^ n \ L J 



I nirx 



>\ r n (23) 



V L 



The boundary term is readily evaluated as A = F J dx'G(x, y; x', 0) = Fy. The lin- 
earization can only be done around the interface of the disorder-free system. This is 
the same as the average interface height of the disordered system only when F is much 
larger than the critical driving force of the underlying pinning-depinning transition. 
For a disorder-free system Eq. JOJ) becomes 

2 J dx'G(x,H Q (ty,x',H (t))d t H Q (t) = FH (t) (24) 

MF 

& d t H (t) = (25) 
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Linearizing Eq. © using H(x,t) = H (t) + h(x,t) leads to 



a^, . 1 , .. „ 1 



-d 2 x h(^ t) - - V (x, H{x, t)) + ~MFh{x, t); (26) 

dx'd y G(x, y\ x', H ) \ y =H h{x,t)d t H ; (27) 

dx'dy<G(x, H ; x', y') \ y ,= Ho h(x', t)d t H ; (28) 

dx'G(x,H ;x',H )d t h(x',t). (29) 



I B + Ic + Id 
Ib 
Ic 
Id 



The derivatives of G, as obtained from the definition of Eq. J23J), are discontinuous at 
y = y' = Ho(t), where the linearization was done. To go around this problem we simply 
set Q(y — y')\ y=y '=Ho(t) — 1/2. We have also performed the same half-strip linearization 
in the case of the spontaneous imbibition, where the half-plane Green's function can 
also be linearized directly [121 121) ■ These two methods yield identical results, i.e. 
linearization and Fourier transformation commute with the half-plane limit. 

The rest of the procedure is then straightforward, and by defining the Fourier series 
representations 



h(p,t) 

v(p,t) 



— I dx cos ) h(x, t); 



L 

pirx 

IT 



r](x,H(x,t)), 



(30) 
(31) 



and projecting Eq. (|26|) to Fourier component p with P p [-] = J dxcos(pnx/L) [•], we 
obtain 



d t H h(p,t) 



1 + e- 2 ^ 



d t h(p,t) 

pix 

2 



1 + e" 2i r 



h(p,t)- V (p,t) + Fh(p,t). 



(32) 



In the Fourier projection of the curvature term, one obtains boundary terms that are 
non-zero at non-zero contact angles, but they are negligible in the limit L — > oo. 
Changing variables to k = pn/L and substituting F = 2d t H we obtain Eq. Q, with 
a discrete wave vector k. Taking the limit L — ► oo while keeping k constant finally 
gives the proper continuum limit. 
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H(x,t) 




(a) 



V> = 



3 y H(y=0) = -F,n(y=0) = 




Figure 1: (a) A schematic geometry and setup of the model. The height of the front 
is described by a single- valued function H(x,t), and the driven boundary condition at 
the reservoir at y = is described by a constant gradient of the chemical potential, (b) 
The profile of the chemical potential (i(x,y) along the y axis at successive time steps 
ti < t 2 ... < t 4 . (c) The profile of the density field 4>(x,y) along the y axis at successive 
time steps corresponding to (b). Note that due to the conservation law these profiles 
have a finite slope in the wet region of the medium, (d) A set of typical rough front 
configurations of a rising interface H(x,t) taken at equal time intervals At = 80. 
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log10(R) Iog10(k) 



Figure 2: (a) Spatial correlation functions G2(r,t) for a system of size L = 256, with 
Aa = 0.3, v = 0.05 at different times. The data are from t = 10 3 to t — 10 4 at 
equal time intervals of t = 2 x 10 3 . In the inset, the crossover lengths £x00 obtained 
from Gi for different velocities are plotted, (b) The structure factor S(k,t) plotted 
against the wave vector k for the same set of parameters as in (a). In the inset, the 
structure factor obtained from the linearized interface equations is plotted against k 
for the corresponding set of parameters, except that L = 512. 



0.5 




Figure 3: Data collapse of the interface width according to the scaling form of Eq. (fTE|) 
for different sets of parameters with L = 256: (i) Aa = 0.3, v = 0.05 and 0.01; (ii) 
Aa = 0.2, v = 0.01,0.005, and 0.002; (iii) Aa = 0.1, v = 0.005,0.002, and 0.001. The 
inset shows the original interface width data. 
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0.4 




Iog10(t) Iog10(t) 

Figure 4: (a) The temporal correlation function C^t) for L = 256, with Aa = 0.2 
and v = 0.015,0.01, and 0.005. In the inset, the corresponding data are shown for the 
linearized interface equations, (b) Temporal correlation functions C q (t) with q — 2,4, 
and 6 for L = 256, with Aa = 0.2 and v = 0.015. In the inset, the the corresponding 
functions are shown for the linearized interface equations. See text for details. 
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